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1.  INTRODUCTION 


An  Important  use  of  time  series  analysis  has  traditionally  been  to 
determine  if  rhythmic  patterns  exist  in  a  time  series  Y.  Udny  Yule  in  1927 
proposed  the  (second  order)  autoregressive  process  as  a  model  for  periodic 
phenomena.  The  purpose  of  this  paper  is  to  investigate  the  autoregressive 
method  for  determining  the  period  of  such  cycles. 

Let  {Y(t),teZ},  Z  the  set  of  integers,  be  a  zero  mean  covariance 
stationary  time  series  with  autocovariance  function  R(v)  =  E(Y(t)Y(t+v) ) , 
vcZ  and  spectral  density  function 

On 

f(u>)  =  £  R(v)e  ivtl1  ,  hte[0,2n]  . 

V— —  ®> 

If  we  assume  that  f  is  bounded  above  and  below  by  positive  constants, 
then  (Masani  (1966))  we  can  write  as  a  limit  in  mean  square 


l  a(j)Y(t-j)  =  e (t)  ,  teZ 
j=0 

where  a(0)  =  1,  the  a(*)  are  real  numbers,  and  (e(t),  teZ)  is  a  white  noise 

series  of  zero  mean  independent,  identically  distributed  random  variables 

2 

with  common  positive  variance  a 

Thus  Y  can  be  written  as  an  infinite  order  autoregressive  process 
(AR(m))  and  f  can  be  written 


c,  \  «  1 

f  ('•>)  =  -y- - - - — 

ztl  I  ,  id).  |  2 

|g(e  )| 

where  the  complex  polynomial  g(z)  =  £ 

unit  circle. 

We  seek  an  estimator  of  a  period 
l.e.  estimate  1  such  that  f * =  0  , 


J 


X 


*=0 


f" 


d)e[0,2ir]  , 

a(j)z^  has  no  zeros  in  or  on  the 


2 it / a)  corresponding  to  a  peak  in  f 

<¥>  -  o  • 
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In  section  2  we  derive  an  estimator  X  of  X  based  on  an  estimated 

AR(p)  approximant  to  f,  i.e.  use  (1.1)  with  g  truncated  at  some  suitable 

<*  A 

degree  p.  Section  3  derives  the  large  sample  distribution  of  X  under  the 

assumption  of  known  finite  autoregressive  order.  The  results  of  a  Monte 
Carlo  simulation  study  illustrating  the  known  order,  unknown  but  finite 
order,  and  infinite  order  cases  are  given  in  section  4  as  well  as  the 
analysis  of  a  real  time -series. 

2.  THE  AUTOREGRESSIVE  METHOD 


Given  a  sample  realization  Y(l),  ....  Y (T)  from  Y,  let  pT(v)  =  R^,(v)/Rt(0), 
v  3  0,  . . . ,  M,  where 

1  T_v 

Rt(v)  =  |  l  Y(t)Y(t+v),  v  =  0,  ....  M  , 

# 

A  A 

and  M  is  an  integer  larger  than  the  approximating  order  p.  Then  p  is  chosen 
as  the  value  of  m  minimizing  Parzen's  (1974)  CAT  criterion 


m  3  0,  M, 


\  1  v  /T- K  —2  .T-m.  — 2 

CAT(m)  3  -  M  t  )oj  "  )0m  * 

j=l  J 

-2 

and  the  a ^  are  found  by  the  Yule-Walker  equations 

l  a  At)  P  T(t-v)  3  S  o  2  ,  v  3  0,  ....  J  <M  , 
1=0  j  T  V  J 


and  6  is  the  Kronecker  delta, 
v 

Thus  the  estimated  AR(p)  spectral  density  f»  is  given  by 


-2 

hM  * 


):  aa(.l)eiJW 
.1-0  P 


,  ta)e[0,2n]  . 


To  estimate  the  period  X  from  f-,  we  need  only  find  the  value  of  w 


maximizing  a  function  of  the  form 
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s(w) 


1 


j=0 


2 


,  uc[0,2tt]  , 


given  p,  a(0)  =  1,  a(l),  a(p)  .  But  this  is  equivalent  to  minimizing 

s  "*■  («o)  =  Hj=0  a(j)eija>)2  which  is  the  spectral  density  of  a  moving  average 
process  of  order  p  (MA(p))  with  parameters  a  and  white  noise  variance  2ir. 
Thus  we  can  write 


s  ^(wj  =  /  Ri(v)e  ^V'  =  Ri(0)  +  2  }  Ri(v)  Cos  vw 

v=-p  v=l 

where  the  inverse  autocovariance  function  (Cleveland  (1972))  Ri(*)  consists 
of  the  autocovariances  corresponding  to  the  MA(p)  process  and  are  given  by 


Ri(v)  =  2ir  y  a(j)a(j+v)  ,  v  =  0,  p  . 

j=0 


Then  given 'an  initial  approximation  w^  to  a  value  w*  of  w  minimizing 
h(w)  =  Ri(v)  Cos  vw  we  use  Newton's  method  to  find  w*  as  the  limit 

of  the  sequence 


h(<iin} 

wn+l  Wn  h'(w  ) 
n 


Ri(v)  Cos  vw 


w  + 
n 


v=l 


vRi(v)  Sin  vw 


v=l 


The  initial  value  w^  is  obtained  in  practice  as  the  frequency  of  a 
relative  maximum  of  f«  evaluated  over  a  suitable  grid  of  frequencies. 

Then  if  K  maxima  of  f»  are  found,  the  estimated  frequencies  are 
labelled  w^,  w^  and  corresponding  periods  by  X.,,  . 

3.  LARGE  SAMPLE  DISTRIBUTION  OF  AN  ESTIMATED  PERIOD 


We  assume  that  Y  is  in  fact  an  AR(p)  process  with  known  p.  Then 


(Parzen  (1961)) 
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/r  (rt  a) 


N  (0, A) 


where  aT  =  (a(l),  a(p)),  aT  =  (a(l),  ....  a(p>),  A  =  has 

(j,k)  element  given  by 
v 


A,.  -  I  [a(j-l-£)a(k-l -l) 
J1C  1=0 

j,k  =  1,  ....  p,  v  =  min  (j-l,k-l), 
of  p  consecutive  Y's. 

From  this  fact  we  have 


-  a ( p+l+£- j ) a ( p+l+£-k) ] 

i 

3 

and  T  is  the  d  x  p  covariance  matrix 
P 


Lemma 

Let  RiT  =(Ri(l),  ....  Ri(p)),  RiT  =  (Ri(l),  Ri(p)),  and  C  be  the 

th  T  1 

p  x  p  matrix  having  v  row  Cv  =  ( Cv  ( 1  >  ,  . ..,  Cv(p))  where  =  Yv+j«(v+j) 

va(j-v),  and  =  1  if  k  =  0,  p  and  0  otherwise. 

Then 

/f(Ri-Ri)  — ^  N(0,CACT)  . 

Proof 


P~v 

/T  (Ri(v)  -  Ri (v) )  =  /r  2tt  £  [«(j)a(j+v)  -  a(j)a(j+v)  ] 

j=0 


P-v 

=  /F  2tt  l  { [a(j)-a(j)  Ja( j+v)4ct (j )  [a(j+v)-a( j+v)  ] } 
j=0 


P~v 

%  /F  l-n  l  { [a(j)-a(j)  ]a(j+v)+a(j)  [a(j+v)-a(j+v)  ] } 
j=0 

P 

l  C  (J)z(J)  , 

J  =  1 

where  z(j)  =  /F  (o(j)-n(D)  ,  j  =  1,  ....  p  and  the  symbol  %  means 
"has  the  same  asymptotic  distribution  as".  Thus  /F  (Ri-Ri^C  /F  (a-a) 
and  the  lemma  follows. 


From  this  we  obtain 
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Theorem 


Let  w  be  the  estimator  obtained  by  the  method  of  section  2  of  the 
true  peak  frequency  w^.  Then 

/t  (w-wQ)  -£-»  N  (0,  bTCACrb) 

where  b(v)  =  v  Sin  vw  /h"(w.)  ,  and  h”(«i)  =  -  j?  v^Ri(v)  Cos  vw  . 

U  v=l 


Proof 


Since  )  vRi(v)  Sin  vw  =  0  and  w  is  a  consistent  estimator  of  wr 


since  f ^  converges  almost  surely  to  we  can  write 


f  v  (Ri(v)-Rl(v) )  Sin  vui 


Q  l  vRi(v)  Sin  vw^  \  vRi(v)[Sin  vw-Sin  vwQ] 
v=l  v=l 


P  .  v— 1 


vRi(v)[Sin  vw-Sin  vw^] 


since 


f  vRi (v) [Sin  vw-Sin  vwQ] 
v=l  _ 


J  vRi(v)[Sin  vw-Sin  vwqI  \  v(Ri  (v)-P.i (v) )  [Sin  vw-Sin  vWq] 
v=l _ .  v=l _ 


f  v(Ri(v)-Ri(v)) (Sin  vw-Sin  VWq) 
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4.  SIMULATION  STUDY  AND  WORKED  EXAMPLE 

To  Illustrate  the  results  of  sections  2  and  3  above  for  known,  unknown 
but  finite,  and  infinite  order  autoregressive  processes  we  performed  a 
simulation  study  in  two  parts:  1)  the  known  and  unknown  order  case,  and 
2)  the  infinite  order  case. 

4.1  Simulation  of  Known  and  Unknown  but  Finite  Order  Case 

We  simulated  on  an  Amdahl  470V/6  computer  at  Texas  A&M  University 

20  series  each  of  length  80,  160,  and  240  for  each  of  the  five  autoregressive 

models  given  in  table  1.  Each  of  the  300  series  was  obtained  by 

Y(t)  =  -][.j|!]a(j)Y(t-j)+e(t),  t  =  p+1,  ...,  T  where  t(l),  ....  e(T)  are 
2 

i id  N(0,a  )  variates  obtained  using  the  Box-Muller  method  (Kennedy  and 
# 

Gentle  (1980),  p.  202)  from  U(0,1)  variates  generated  by  a  composite 

generator  (Kennedy  and  Gentle  (1980),  p.  162),  and  starting  values  Y(l), 

...,  Y(p)  were  obtained  as  a  N  (0,T  )  random  variate  where  F  is  obtained 

p  ■  p  P 

by  the  algorithm  of  Pagano  (1973). 

Estimators  of  the  spectral  density  for  each  series  were  determined 
in  four  ways: 

1)  p,YW:  the  Yule-Walker  equations  were  solved  via  Levinson's 
algorithm  (Levinson  (1947))  to  determine  estimates  of  the  parameters 
of  the  AR(p)  process. 

2 

2)  p,LS:  least  squares  estimators  of  u(l),  ...,  a(p),  a  were  obtained 
via  a  Gram-Sehmidt  orthogonal ization  procedure  applied  to  the 
regression  model  Y  =  -Xn+c  where  Y^  =  (Y(p+1),  ...,  Y(T)  , 

aT  =  (a ( 1) ,  ...,  a(p) ) ,  X.^  =  Y(p+j-k),  j  =  1,  ...,  T-p  ,  k*l, 

...,  p  ,  and  r. F  =  (f(p+l),  ....  f(T))  . 
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3)  p,YW:  The  CAT  criterion  determined  order  (with  6^  obtained  via 
Levinson's  algorithm)  parameters  estimated  by  the  Yule-Walker  equations. 

4)  p,LS:  Same  as  p,YW  except  parameters  estimated  by  least  squares. 

Figure  A  contains  plots  of  the  true  log  spectra  for  each  of  the  five 

models,  while  table  1  lists  the  true  periods  (obtained  from  the  coefficients 

2 

by  the  algorithm  of  section  2)  for  each  model  as  well  as  Ta»  .  In  each 

case  the  peak  of  interest  is  the  one  for  smallest  frequency.  We  note  that 

models  1,  2,  and  5  are  very  similar  except  for  order  and  location  of  peak 

2 

while  model  4  has  a  single  broad  peak  (reflected  in  the  size  of  To^). 

Model  3  is  of  particular  interest  since  it  appears  that  h"(8.951)  %  0 

2 

(note  how  this  is  reflected  in  To^  ). 

Thus  it  appears  that  these  models  are  representative  of  a  wide  class 
of  models  and  that  one  pathological  case  (model  3)  has  been  included. 

The  purpose  of  this  part  of  the  simulation  study  is  twofold: 

1)  In  the  p,YW  and  p,LS  case  to  see  if  the  methodology  of  section  2 
agrees  with  the  theory  of  section  3  and  to  compare  the  performance  of  YW 
and  LS,  particularly  on  model  3. 

2)  To  investigate  the  distribution  of  X  in  the  finite  but  unknown 
order  case. 


Table  2  contains  results  for  the  p,YW  and  p,LS  cases  while  figures 


B-F  display  the  20  p,YW  and  p,I.S  estimated  spectra  for  the  three  sample 
sizes  for  each  of  the  five  models.  A  comparison  of  the  spectra  for  the 
YW  and  LS  estimators  shows  they  are  almost  indistinguishable  except  that 
LS  has  a  wider  spectral  range  in  models  2,  3,  and  5  and  that  the  peak  in 
model  3  is  estimated  somewhat  differently  by  the  two  methods.  This  is 
reflected  in  table  2  where  it  is  seen  that  LS  decomposes  this  peak  into 
two  peaks  in  4,  4,  and  3  series  (T  =  80,  T  ■  160,  T  *  240)  while  YW  does 
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tliis  ill  3,  ],  and  1  series.  Tims  except  for  this  pathological  model  we 
find  no  large  difference  in  the  two  estimation  procedures.  Also  table  2 
shows  the  expected  adherence  of  the  simulated  data  to  the  theory  of 
section  3. 

Table  3  and  figures  G-K  report  the  results  of  the  p,YW  analysis  of 
the  300  simulated  autoregressive  processes.  The  results  indicate  that 
the  method  estimates  periods  remarkably  well.  In  fact,  even  though  45% 
of  the  series  had  the  wrong  order  determined,  the  results  of  table  3  are 
remarkably  similar  to  those  of  table  2. 

4.2  Simulation  of  Infinite  Order  Case 

Twenty  series  of  lengths  80,  160,  and  240  were  generated  from  the 
MA ( 8 )  model  Y(t)  =  e(t)  +  .309e(t-l)  -  ,0748t(t-2  -  ,0113e(t-3)  -  .0853e(t-4) 
-  ,0552r. (t-5)  +  .0084e(t-6)  +  .4621e(t-7)  +  .288c(t-8)  .  Figure  L  contains 
the  true  MA(8)  log  spectrum  and  the  estimated  autoregressive  spectra  for 
each  of  the  60  simulated  series.  Table  4  describes  the  analysis  of  each  of 
the  60  series. 

Inspection  of  table  4  raises  two  points: 

1)  The  AR  orders  chosen  appear  to  increase  with  sample  size. 

2)  This  increase  leads  to  estimated  models  which  make  two  (or  even 
three  in  one  case)  peaks  out  of  the  single  broad  peak. 

Thus  it  appears  that  if  the  process  Y  cannot  be  adequately  represented 
as  a  finite  order  AK  process  (in  particular  if  the  process  has  very  broad 
peaks),  the  proposed  method  may  not  be  satisfactory. 

We  note  however  that  if  instead  of  using  approximating  AR  schemes 
one  used  approximating  MA  schemes,  this  would  make  virtually  no  change  in 
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Lho  algorithm  and  only  the  maLrlx  A  in  the  asymptotic  covariance  of  X 
would  need  to  be  changed  in  theorem  1  to  the  asymptotic  covariance  matrix 
of  maximum  likelihood  estimators  of  MA  parameter  estimators.  We  note 
however  that  these  methods  cannot  be  extended  to  the  mixed  AR-MA  processes 
since  neither  their  spectra  nor  its  reciprocal  can  be  written  as  a  finite 
degree  trigonometric  polynomial.  In  this  paper  we  have  confined  our 
attention  to  approximating  AR  schemes  because  of  their  computational  at¬ 
tractiveness  and  widespread  use. 

4.3  Analysis  of  Hormone  Levels  Data 

We  consider  the  levels  of  luteinizing  hormone  (LH)  in  a  cow  as 
measured  at  10  minute  intervals  for  a  24  hour  period  (Rahe  et  al  (1980)). 

A  plot  of  (he  data  is  given  in  figure  M.  Figure  N  displays  the  log 
spectra  of  the  p,YW  determined  AR(13)  model 

Y(t)  -  .2736  Y(t-l)  +  .0769  Y(t-2)  +  .0046  Y(t-3)  +  .1576  Y(t-4) 

-  .0131  Y ( t— 5 )  +  .1486  Y(t-6)  -  .4795  Y(t-7)  -  .1942  Y(t-8) 

+  .0546  Y ( t-9)  +  .0876  Y(t-10)  -  .1546  Y(t-ll  +  .0756  Y(t-12) 

-  .1926  Y ( t-13)  =  e(t)  . 

This  model  gives  X  =  7.278  time  intervals  with  an  estimated  standard 
error  s<^  =  .116.  Thus  a  large  sample  confidence  interval  for  X  is  given 
by  (70.46  min,  75.10  min)  . 
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1.  Autoregressive  Models  used  in  Simulation  Study 


Model 

Number 

P 

1 

Tos 

Coefficients 

1 

7 

5.900 

3.120 

-.6381 

-.0140 

.0295 

.0272 

-.1587 

-.4581 

.3372 

2 

7 

7.398 

4.240 

-.2004 

.1158 

.1679 

.1056 

.1357 

.0566 

-.3899 

3 

9 

8.951 

924.0 

-.4335 

.3740 

.1071 

.3213 

.2071 

.2089 

.1942 

-.2828 

.3656 

4 

4 

• 

6.900 

13.44 

-.1477 

.0601 

.3565 

.1964 

5 

10 

7.586 

4.80 

-.3190 

-.0175 

.1803 

.0340 

.1832 

0423 


1393 


2463 


1365 


1974 


2.  Simulation  Results  for  Known  Order  Case4 


Model  1 

Model  2 

Model  3 

Model  4 

Model  5 

p 

7 

7 

9 

4 

10 

X 

5.900 

7.398 

8.951 

6.900 

7.586 

°£.«o 

.039 

.053 

11.550 

.168 

.036 

°Ll60 

.020 

.026 

5.775 

.084 

.018 

°X.240 

.013 

.018 

3.850 

.056 

.012 

XAVE 

5.912,5.937 

7.410,7.412 

8.776,8.82 

6.897,6.891 

7.480,7.504 

2 

SX 

.056,  .041 

.198,  .173 

.624,  .736 

.158,  .148 

.077,  .071 

N 

00 

o 

p 

CIC 

.950,  .900 

.750,  .  750 

1.000,1.000 

1.000,  .950 

.850,  .800 

NBC 

• 

3  ,4 

'^AVE 

5.983,5.914 

7.474,7.474 

8.997.8.961 

6.901,6.906 

7.598,7.606 

2 

3X 

>160 

.041,  .019 

.034,  .033 

.500,  .504 

.052,  .054 

.030,  .028 

CIC 

1.000,  .950 

.900,  .900 

1.000,1.000 

.950,  .950 

.850,  .950 

NB 

1  .4 

*AVE 

5.871,5.874 

7.431,7.428 

8.889,8.82 

6.832,6.832 

7.556,7.564 

2 

.023,  .021 

.021,  .021 

.575,  .652 

.042,  .041 

.009,  .009 

-240 

CIC 

.900,  .900 

1.000,1.000 

1.000,1.000 

1.000,1.000 

1.000,1.000 

NB 

1  .3 

3  First  number  In  a  pair  Is  for  Yule-Walker  estimator,  second  is  for  Least 
Squares  estimator. 

b  Confidence  Interval  Coverage.  (Proportion  of  the  20  X  within  2^  of  X)  . 

C  Number  of  Series  producing  dual  peaks  In  model  3.  These  series  not  Included 
In  calculation  of  £  ,  sj:  ,  CIC. 
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3.  Simulation  Results  for  Unknown  Order  Case*1 


Model  1 

Model  2 

Model  3 

Model  4 

Model  5 

P 

7 

7 

9 

4 

10 

A 

5.900 

7.398 

8.951 

6.900 

7.586 

°X,80 

.039 

.053 

11.550 

.168 

.036 

aS,160 

.020 

.026 

5.775 

.084 

1 

.018 

aS,240 

*  .013 

.018 

3.850 

.056 

.012 

P  AVE  * 9  p 

7.050,  .366 

6.650,9.621 

7.650,5.713 

3.850,1.082 

7.150,10.239 

T-80 

pmin,pmax,pccC 

6,8, .65 

0,13, .40 

3,13, .43 

3, 7, .50 

0,11, .15 

X,s“- 

5.944,  .061 

7.384^  .310 

8.634,  .832 

6.786,  .534 

7.578fa,  .151 

CIC,  NB 

.95 

1.C0 

1.00  ,3 

.85 

.80 

PAVE ’ 8p  ' 

7.450,  .471 

7.950,3.418 

9.600,  .568 

5.050,2.050 

9.800,4.695 

PMIN*PMAX,PCC 

7,9, .65 

7, 13, .70 

9, 11, .55 

3,8, .35 

3,14, .60 

T-160 

5.919,  .016 

7.479, .049 

8.923,  .557 

6.811,  .395 

7.632,  .057 

CIC, NB 

1.00 

.85 

1.00  ,3 

.65 

.85 

2 

PAVE,8p 

8.000,2.947 

7.550,1.520 

10.30,3.800 

4.900,2.305 

10.45,  .787 

PMIN*PMAX,PCC 

7,13, .60 

7, 11, .80 

9,15,-50 

3.8, .60 

10, 13, .75 

T-240 

A.sJ 

5.878,  .045 

7.449,  .026 

9.055,  .592 

6.830,  .010 

7.550,  .012 

CIC.NB 

.85 

.95 

1.00  ,3 

.80 

.95 

a  Yule-Walker  estimators  used. 

fa  Series  having  no  peaks  not  Included  In  calculations. 
C.  ProDortlon  choosinR  correct  order. 
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4.  Simulation  Results  for  MA(8)  Model  4 


T  - 

80 

T  -  160 

T  -  240 

Series 

P 

Periods 

P 

Periods 

P 

Periods 

Number 

1 

14  9. 

204,5.888 

14 

8.552,5.912 

15 

8.159 

2 

7 

7.141 

7 

7.578 

9 

7.473 

3 

8 

6.407 

13 

8.546,5.815 

19 

7.736,5.869 

4 

7 

7.881 

12 

8.470,6.293 

7 

7.409 

5 

0 

— 

7 

7.964 

7 

7.S29 

6 

8 

7.238 

14 

9.357,5.997 

9 

7.540 

7 

2 

5.382 

7 

7.559 

14 

8.871,6.278 

8 

7 

7.452 

7 

7.324 

7 

7.351 

9 

e 

7 

8.015 

7 

7.110 

14 

9.117,6.351 

10 

7 

7.312 

7 

7.817 

8 

7.209 

11 

7 

7.386 

7 

7.097 

10 

7.896 

12 

4 

6.328 

7 

7.541 

7 

7.585 

13 

7 

7.678 

7 

7.133 

7 

7.111 

14 

7 

7.508 

7 

7.408 

14 

8.824,5.906 

15 

7 

6.916 

8 

6.589 

15 

7.563,5.807 

16 

7 

7.216 

14 

8.500,6.280 

14 

7.788 

17 

2 

5.021 

7 

7.750 

18 

9.785,7.393,5.529 

18 

8 

6.847 

7 

7.255 

15 

9.181,6.367 

19 

7 

7.365 

7 

7.060 

9 

8.405 

20 

7 

6.787 

7 

7.335 

14 

8.937,6.019 

Mea«k 

6.S00 

6.993 

8.650 

7.368 

11 

.600 

7.621 

iriance 

8.158 

.591 

8.134 

.115 

15 

.832 

.133 

a  YW  estimators  used,  true  period  Is  7.334. 

6  These  calculations  for  periods  do  not  include  peeks  or  multi-peeked  series. 
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